Data collection was divided into Canopy and Understory

Dictionary of terms:

plant_part = canopy or understory (under),

nitrate = nutrient dose, which includes phosphate,

d9_growth_percent = calculated growth from final weight-initial weight divided by initial,

mean_mins28 = mean of summed time in minutes for 9 day period each bin in canopy or understory experienced temperatures 28C or higher

Midway Experiment
Midway Experiment

Prepare the data

Load packages

library(lme4)
library(lmerTest)
library(effects)
library(car)
library(MuMIn)
library (dplyr)
library(emmeans)
library(DHARMa)
library(performance)
library(patchwork)
library(rstatix)
#for plots and tables
library(ggplot2)
library(ggpubr)
library(forcats)
library(RColorBrewer)
library(tidyverse)
library(sjPlot)
library(sjmisc)
library(gt)
library(purrr)
library(stringr)
library(tidyr)
library(piecewiseSEM)
library(easystats)
library(magrittr)

Open weight dataset

mid_growth <- read.csv("../../../data/midway_2024/input/midway_growth.csv")

Calculate growth rate from initial and final weights, make new columns

mid_growth$d9_growth_percent <- round(((mid_growth$final_weight - mid_growth$initial_weight) / mid_growth$initial_weight * 100), digits = 2)
mid_growth$d5_growth_percent <- round(((mid_growth$d5_weight - mid_growth$initial_weight) / mid_growth$initial_weight * 100), digits = 2)

Add a new column that gets rid of characters in ID column

mid_growth$plant_ID <- as.factor(substr(mid_growth$ID, 4, 5))

Make a new column for run using the third value in the ID, then change the 1 to 2 and 2 to 3

mid_growth$run <- substr(mid_growth$ID, 3, 3)
mid_growth$run <- as.factor(ifelse(mid_growth$run == 1, 2, 3)) #actual run numbers are 2 and 3

Assign treatment as characters from integers then to factors

mid_growth$nitrate <- as.factor(as.character(mid_growth$treatment))

Get rid of the ppt in salinity

mid_growth$salinity <- substr(mid_growth$salinity, 1, 2)

Combine nitrate and salinity for a treatment number

mid_growth <- mid_growth %>%
  mutate(treatment = case_when(
    nitrate == 1 & salinity == 35 ~ "1",
    nitrate == 1 & salinity == 28 ~ "2",
    nitrate == 2 & salinity == 35 ~ "3",
    nitrate == 2 & salinity == 28 ~ "4",
    nitrate == 3 & salinity == 35 ~ "5",
    nitrate == 3 & salinity == 28 ~ "6",
    nitrate == 4 & salinity == 35 ~ "7",
    nitrate == 4 & salinity == 28 ~ "8",
  ))

Create subsets for the plots

Use Day 9 for final analysis

canopy_g <- subset(mid_growth, plant_part == "canopy" & d9_growth_percent != 42.45)
canopy_g$treatment_graph[canopy_g$treatment == 1] <- "1) 0.5umol/35 ppt"
canopy_g$treatment_graph[canopy_g$treatment == 2] <- "2) 0.5umol/28 ppt"
canopy_g$treatment_graph[canopy_g$treatment == 3] <- "3) 2umol/35 ppt" 
canopy_g$treatment_graph[canopy_g$treatment == 4] <- "4) 2umol/28 ppt"
canopy_g$treatment_graph[canopy_g$treatment == 5] <- "5) 4umol/35 ppt"
canopy_g$treatment_graph[canopy_g$treatment == 6] <- "6) 4umol/28 ppt"
canopy_g$treatment_graph[canopy_g$treatment == 7] <- "7) 8umol/35 ppt"
canopy_g$treatment_graph[canopy_g$treatment == 8] <- "8) 8umol/28 ppt"
under_g <- subset(mid_growth, plant_part == "under" & d9_growth_percent != -62.82 & 
                    d9_growth_percent != -32.50) #removing one influential datapoint
under_g$treatment_graph[under_g$treatment == 1] <- "1) 0.5umol/35 ppt"
under_g$treatment_graph[under_g$treatment == 2] <- "2) 0.5umol/28 ppt"
under_g$treatment_graph[under_g$treatment == 3] <- "3) 2umol/35 ppt" 
under_g$treatment_graph[under_g$treatment == 4] <- "4) 2umol/28 ppt"
under_g$treatment_graph[under_g$treatment == 5] <- "5) 4umol/35 ppt"
under_g$treatment_graph[under_g$treatment == 6] <- "6) 4umol/28 ppt"
under_g$treatment_graph[under_g$treatment == 7] <- "7) 8umol/35 ppt"
under_g$treatment_graph[under_g$treatment == 8] <- "8) 8umol/28 ppt"

Add new column to subsets from “time over 28” summary dataset

time_over_28 <- read_csv("../../../data/midway_2024/transformed/mins_28_plant_part.csv") #load dataset

Keep only relevant columns, convert to factors before combining datasets so data types match

mean_mins28 <- time_over_28 %>%
  select(plant_part, run, mean_mins28_plant_part) %>% #keep only relevant columns of data
  mutate(run = as.factor(run), mean_mins28_plant_part = as.factor(mean_mins28_plant_part))

Join the time over 28 dataset to the growth dataset using run and plant_part as matches CANOPY

canopy_g <- canopy_g %>%
  left_join(mean_mins28, by = c("run", "plant_part")) #join the datasets

canopy_g <- canopy_g %>%
  rename(mean_mins28 = mean_mins28_plant_part) #name is too long
glimpse(canopy_g)
## Rows: 96
## Columns: 14
## $ ID                <chr> "ct101", "ct102", "ct103", "ct104", "ct105", "ct106"…
## $ treatment         <chr> "1", "2", "3", "4", "5", "6", "7", "8", "1", "2", "3…
## $ salinity          <chr> "35", "28", "35", "28", "35", "28", "35", "28", "35"…
## $ plant_part        <chr> "canopy", "canopy", "canopy", "canopy", "canopy", "c…
## $ initial_weight    <dbl> 0.3964, 0.3709, 0.4576, 0.4222, 0.3401, 0.3468, 0.42…
## $ d5_weight         <dbl> 0.4140, 0.3949, 0.4280, 0.4321, 0.3892, 0.3963, 0.43…
## $ final_weight      <dbl> 0.4159, 0.4040, 0.4143, 0.4307, 0.4150, 0.4028, 0.47…
## $ d9_growth_percent <dbl> 4.92, 8.92, -9.46, 2.01, 22.02, 16.15, 12.14, -42.45…
## $ d5_growth_percent <dbl> 4.44, 6.47, -6.47, 2.34, 14.44, 14.27, 2.29, -10.94,…
## $ plant_ID          <fct> 01, 02, 03, 04, 05, 06, 07, 08, 09, 10, 11, 12, 13, …
## $ run               <fct> 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2…
## $ nitrate           <fct> 1, 1, 2, 2, 3, 3, 4, 4, 1, 1, 2, 2, 3, 3, 4, 4, 1, 1…
## $ treatment_graph   <chr> "1) 0.5umol/35 ppt", "2) 0.5umol/28 ppt", "3) 2umol/…
## $ mean_mins28       <fct> 2825, 2825, 2825, 2825, 2825, 2825, 2825, 2825, 2825…

Join the time over 28 dataset to the growth dataset using run and plant_part as matches UNDERSTORY

under_g <- under_g %>%
  left_join(mean_mins28, by = c("run", "plant_part")) #join the datasets

under_g <- under_g %>%
  rename(mean_mins28 = mean_mins28_plant_part) #name is too long
glimpse(under_g)
## Rows: 94
## Columns: 14
## $ ID                <chr> "ct149", "ct150", "ct151", "ct152", "ct153", "ct154"…
## $ treatment         <chr> "1", "2", "3", "4", "5", "6", "7", "8", "1", "2", "3…
## $ salinity          <chr> "35", "28", "35", "28", "35", "28", "35", "28", "35"…
## $ plant_part        <chr> "under", "under", "under", "under", "under", "under"…
## $ initial_weight    <dbl> 0.4615, 0.4468, 0.3672, 0.4865, 0.4462, 0.4084, 0.48…
## $ d5_weight         <dbl> 0.4983, 0.4398, 0.3940, 0.4496, 0.4498, 0.4422, 0.46…
## $ final_weight      <dbl> 0.5332, 0.4092, 0.3755, 0.4294, 0.4265, 0.4649, 0.40…
## $ d9_growth_percent <dbl> 15.54, -8.42, 2.26, -11.74, -4.42, 13.83, -16.43, 0.…
## $ d5_growth_percent <dbl> 7.97, -1.57, 7.30, -7.58, 0.81, 8.28, -3.42, 2.08, 0…
## $ plant_ID          <fct> 49, 50, 51, 52, 53, 54, 55, 56, 57, 58, 59, 61, 62, …
## $ run               <fct> 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2…
## $ nitrate           <fct> 1, 1, 2, 2, 3, 3, 4, 4, 1, 1, 2, 3, 3, 4, 4, 1, 1, 2…
## $ treatment_graph   <chr> "1) 0.5umol/35 ppt", "2) 0.5umol/28 ppt", "3) 2umol/…
## $ mean_mins28       <fct> 1398, 1398, 1398, 1398, 1398, 1398, 1398, 1398, 1398…

Analysis – Canopy

Run linear mixed-effects model and run checks for fit

Make histograms

canopy_g %>% 
ggplot(aes(d9_growth_percent)) + #day 9
  geom_histogram(binwidth=5, fill = "goldenrod1", color = "black", linewidth = 0.25, alpha = 0.85) +
  theme_bw()

canopy_g %>% ggplot(aes(d5_growth_percent)) + #day 5
  geom_histogram(binwidth=5, fill = "goldenrod1", color = "black", linewidth = 0.25, alpha = 0.85) +
  theme_bw()

Run the model

Mean_min28 used as random effect in lieu of run, but neither had an effect, also rlc_order had no effect

growth_model_canopy <- lmer(formula = d9_growth_percent ~ salinity + nitrate +
                            (1 | plant_ID), data = canopy_g, REML = TRUE)

Check model fit

isSingular(growth_model_canopy)
## [1] FALSE
hist(resid(growth_model_canopy))

plot(resid(growth_model_canopy) ~ fitted(growth_model_canopy))

qqnorm(resid(growth_model_canopy))
qqline(resid(growth_model_canopy))

Check the performance of the model for dataset: canopy

performance::check_model(growth_model_canopy)

r.squaredGLMM(growth_model_canopy)
##             R2m       R2c
## [1,] 0.09262105 0.3014886
summary(growth_model_canopy)
## Linear mixed model fit by REML. t-tests use Satterthwaite's method [
## lmerModLmerTest]
## Formula: d9_growth_percent ~ salinity + nitrate + (1 | plant_ID)
##    Data: canopy_g
## 
## REML criterion at convergence: 738.7
## 
## Scaled residuals: 
##      Min       1Q   Median       3Q      Max 
## -3.03933 -0.63855  0.03797  0.56284  2.36387 
## 
## Random effects:
##  Groups   Name        Variance Std.Dev.
##  plant_ID (Intercept)  39.51    6.286  
##  Residual             132.12   11.495  
## Number of obs: 96, groups:  plant_ID, 48
## 
## Fixed effects:
##             Estimate Std. Error     df t value Pr(>|t|)  
## (Intercept)    5.639      3.316 43.000   1.701   0.0963 .
## salinity35     5.960      2.966 43.000   2.009   0.0508 .
## nitrate2      -3.383      4.195 43.000  -0.807   0.4243  
## nitrate3       2.340      4.195 43.000   0.558   0.5798  
## nitrate4      -5.145      4.195 43.000  -1.227   0.2266  
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Correlation of Fixed Effects:
##            (Intr) slnt35 nitrt2 nitrt3
## salinity35 -0.447                     
## nitrate2   -0.632  0.000              
## nitrate3   -0.632  0.000  0.500       
## nitrate4   -0.632  0.000  0.500  0.500

View random effects levels

tab_model(growth_model_canopy, show.intercept = TRUE, show.se = TRUE, show.stat = TRUE, show.df = TRUE, show.zeroinf = TRUE)
  d9_growth_percent
Predictors Estimates std. Error CI Statistic p df
(Intercept) 5.64 3.32 -0.95 – 12.23 1.70 0.093 89.00
salinity [35] 5.96 2.97 0.07 – 11.85 2.01 0.048 89.00
nitrate [2] -3.38 4.19 -11.72 – 4.95 -0.81 0.422 89.00
nitrate [3] 2.34 4.19 -5.99 – 10.68 0.56 0.578 89.00
nitrate [4] -5.15 4.19 -13.48 – 3.19 -1.23 0.223 89.00
Random Effects
σ2 132.12
τ00 plant_ID 39.51
ICC 0.23
N plant_ID 48
Observations 96
Marginal R2 / Conditional R2 0.093 / 0.301
plot(allEffects(growth_model_canopy))

Construct null model to perform likelihood ratio test (REML must be FALSE)

canopy_nitrate_null <- lmer(formula = d9_growth_percent ~ salinity + 
                              (1 | plant_ID), data = canopy_g, REML = FALSE)
canopy_model2 <- lmer(formula = d9_growth_percent ~ nitrate + salinity +
                        (1 | plant_ID), data = canopy_g, REML = FALSE)
anova(canopy_nitrate_null, canopy_model2)
## Data: canopy_g
## Models:
## canopy_nitrate_null: d9_growth_percent ~ salinity + (1 | plant_ID)
## canopy_model2: d9_growth_percent ~ nitrate + salinity + (1 | plant_ID)
##                     npar    AIC    BIC  logLik deviance  Chisq Df Pr(>Chisq)
## canopy_nitrate_null    4 770.61 780.86 -381.30   762.61                     
## canopy_model2          7 772.50 790.45 -379.25   758.50 4.1109  3     0.2497
canopy_salinity_null <- lmer(formula = d9_growth_percent ~ nitrate + 
                               (1 | plant_ID), data = canopy_g, REML = FALSE)
canopy_model3 <- lmer(formula = d9_growth_percent ~ nitrate + salinity + 
                        (1 | plant_ID), data = canopy_g, REML = FALSE)
anova(canopy_salinity_null, canopy_model3)
## Data: canopy_g
## Models:
## canopy_salinity_null: d9_growth_percent ~ nitrate + (1 | plant_ID)
## canopy_model3: d9_growth_percent ~ nitrate + salinity + (1 | plant_ID)
##                      npar   AIC    BIC  logLik deviance Chisq Df Pr(>Chisq)  
## canopy_salinity_null    6 774.8 790.19 -381.40    762.8                      
## canopy_model3           7 772.5 790.45 -379.25    758.5 4.308  1    0.03793 *
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Plot canopy growth

canopy_growth_plot <- canopy_g %>% 
  ggplot(aes(treatment_graph, d9_growth_percent)) + 
  geom_boxplot(size=0.5) + 
  geom_point(alpha = 0.75, size = 5, aes(color = salinity), position = "jitter", show.legend = TRUE) + 
  labs(x="Nitrate", y= "9-Day Growth (%)", title= "A", subtitle = "Chondria tumulosa -- Canopy") + 
  scale_x_discrete(labels = c("0.5 μmol N", "0.5 μmol N", "2 μmol N", "2 μmol N", 
                              "4 μmol N", "4 μmol N", "8 μmol N", "8 μmol N")) + 
  ylim(-45, 50) + stat_mean() + 
  scale_color_manual(values = c("goldenrod1", "darkgoldenrod3")) +
  geom_hline(yintercept=0, color = "red", linewidth = 0.5, alpha = 0.5) +
  theme_bw() +
  theme(legend.position.inside = c(0.90,0.90), plot.title = element_text(face = "bold", vjust = -15, hjust = 0.05), 
        plot.subtitle = element_text(face = "italic", size = 14, vjust = -20, hjust = 0.05))
canopy_growth_plot

Analysis – Understory

Make a histogram of the growth rate data

under_g %>% ggplot(aes(d9_growth_percent)) + #day 9
  geom_histogram(binwidth=5, fill = "#990066", color = "black", linewidth = 0.25, alpha = 0.85) +
  theme_bw()

under_g %>% ggplot(aes(d5_growth_percent)) + #day 5
  geom_histogram(binwidth=5, fill = "#990066", color = "black", linewidth = 0.25, alpha = 0.85) +
  theme_bw()

Mean_min28 used as random effect in lieu of run, RLC order had no effect

growth_model_under <- lmer(formula = d9_growth_percent ~ salinity + nitrate +
                              (1 | plant_ID) + (1 | mean_mins28), data = under_g, REML = TRUE)

isSingular(growth_model_under)
## [1] FALSE
hist(resid(growth_model_under))

plot(resid(growth_model_under) ~ fitted(growth_model_under))

qqnorm(resid(growth_model_under))
qqline(resid(growth_model_under))

Check the performance of the model for dataset: under

performance::check_model(growth_model_under)

r.squaredGLMM(growth_model_under)
##             R2m       R2c
## [1,] 0.07896071 0.2748286
summary(growth_model_under)
## Linear mixed model fit by REML. t-tests use Satterthwaite's method [
## lmerModLmerTest]
## Formula: d9_growth_percent ~ salinity + nitrate + (1 | plant_ID) + (1 |  
##     mean_mins28)
##    Data: under_g
## 
## REML criterion at convergence: 653.7
## 
## Scaled residuals: 
##      Min       1Q   Median       3Q      Max 
## -2.50427 -0.51736 -0.02634  0.54235  2.21577 
## 
## Random effects:
##  Groups      Name        Variance Std.Dev.
##  plant_ID    (Intercept)  0.9344  0.9666  
##  mean_mins28 (Intercept) 18.7592  4.3312  
##  Residual                72.9123  8.5389  
## Number of obs: 94, groups:  plant_ID, 48; mean_mins28, 2
## 
## Fixed effects:
##             Estimate Std. Error     df t value Pr(>|t|)   
## (Intercept)   13.207      3.646  1.704   3.623  0.08643 . 
## salinity35    -1.844      1.785 42.742  -1.033  0.30730   
## nitrate2      -7.434      2.524 42.720  -2.946  0.00519 **
## nitrate3      -3.589      2.524 42.720  -1.422  0.16218   
## nitrate4      -3.287      2.496 41.754  -1.317  0.19510   
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Correlation of Fixed Effects:
##            (Intr) slnt35 nitrt2 nitrt3
## salinity35 -0.245                     
## nitrate2   -0.335 -0.015              
## nitrate3   -0.335 -0.015  0.490       
## nitrate4   -0.342  0.000  0.495  0.495

View random effects levels

tab_model(growth_model_under, show.intercept = TRUE, show.se = TRUE, show.stat = TRUE, show.df = TRUE, show.zeroinf = TRUE)
  d9_growth_percent
Predictors Estimates std. Error CI Statistic p df
(Intercept) 13.21 3.65 5.96 – 20.46 3.62 <0.001 86.00
salinity [35] -1.84 1.78 -5.39 – 1.70 -1.03 0.304 86.00
nitrate [2] -7.43 2.52 -12.45 – -2.42 -2.95 0.004 86.00
nitrate [3] -3.59 2.52 -8.61 – 1.43 -1.42 0.159 86.00
nitrate [4] -3.29 2.50 -8.25 – 1.68 -1.32 0.191 86.00
Random Effects
σ2 72.91
τ00 plant_ID 0.93
τ00 mean_mins28 18.76
ICC 0.21
N plant_ID 48
N mean_mins28 2
Observations 94
Marginal R2 / Conditional R2 0.079 / 0.275
plot(allEffects(growth_model_under))

Construct null model to perform likelihood ratio test REML must be FALSE

under_nitrate_null <- lmer(formula = d9_growth_percent ~ salinity + 
                              (1 | plant_ID) + (1 | mean_mins28), data = under_g, REML = FALSE)
under_model2 <- lmer(formula = d9_growth_percent ~ nitrate + salinity +
                        (1 | plant_ID) + (1 | mean_mins28), data = under_g, REML = FALSE)
anova(under_nitrate_null, under_model2)
## Data: under_g
## Models:
## under_nitrate_null: d9_growth_percent ~ salinity + (1 | plant_ID) + (1 | mean_mins28)
## under_model2: d9_growth_percent ~ nitrate + salinity + (1 | plant_ID) + (1 | mean_mins28)
##                    npar    AIC    BIC  logLik deviance Chisq Df Pr(>Chisq)  
## under_nitrate_null    5 689.38 702.10 -339.69   679.38                      
## under_model2          8 686.74 707.08 -335.37   670.74 8.644  3    0.03442 *
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
under_salinity_null <- lmer(formula = d9_growth_percent ~ nitrate + (1 | plant_ID) + (1 | mean_mins28), data = under_g, REML = FALSE)
under_model3 <- lmer(formula = d9_growth_percent ~ nitrate + salinity + (1 | plant_ID) + (1 | mean_mins28), data = under_g, REML = FALSE)
anova(under_salinity_null, under_model3)
## Data: under_g
## Models:
## under_salinity_null: d9_growth_percent ~ nitrate + (1 | plant_ID) + (1 | mean_mins28)
## under_model3: d9_growth_percent ~ nitrate + salinity + (1 | plant_ID) + (1 | mean_mins28)
##                     npar    AIC    BIC  logLik deviance  Chisq Df Pr(>Chisq)
## under_salinity_null    7 685.88 703.68 -335.94   671.88                     
## under_model3           8 686.74 707.08 -335.37   670.74 1.1413  1     0.2854

Plot understory growth

under_growth_plot <- under_g %>% 
  ggplot(aes(treatment_graph, d9_growth_percent)) + 
  geom_boxplot(size=0.5) + 
  geom_point(alpha = 0.75, size = 5, aes(color = salinity), position = "jitter", show.legend = TRUE) + 
  labs(x="Nitrate", y= "9-Day Growth (%)", title= "B", subtitle = "Chondria tumulosa -- Understory") + 
  scale_x_discrete(labels = c("0.5 μmol N", "0.5 μmol N", "2 μmol N", "2 μmol N", 
                              "4 μmol N", "4 μmol N", "8 μmol N", "8 μmol N")) + 
  ylim(-35, 40) + stat_mean() + 
  scale_color_manual(values = c("maroon", "maroon2")) +
  geom_hline(yintercept=0, color = "red", linewidth = 0.5, alpha = 0.5) +
  theme_bw() +
  theme(legend.position.inside = c(0.90,0.90), plot.title = element_text(face = "bold", vjust = -15, hjust = 0.05), 
        plot.subtitle = element_text(face = "italic", size = 14, vjust = -20, hjust = 0.05))
under_growth_plot

Summarize the means for nitrate and salinity

mid_growth %>% group_by(nitrate, plant_part) %>% summarise_at(vars(d9_growth_percent), list(mean = mean))
## # A tibble: 8 × 3
## # Groups:   nitrate [4]
##   nitrate plant_part  mean
##   <fct>   <chr>      <dbl>
## 1 1       canopy      8.62
## 2 1       under      12.3 
## 3 2       canopy      5.24
## 4 2       under       2.12
## 5 3       canopy     11.0 
## 6 3       under       7.07
## 7 4       canopy      3.47
## 8 4       under       9.00
mid_growth %>% group_by(salinity, plant_part) %>% summarise_at(vars(d9_growth_percent), list(mean = mean))
## # A tibble: 4 × 3
## # Groups:   salinity [2]
##   salinity plant_part  mean
##   <chr>    <chr>      <dbl>
## 1 28       canopy      4.09
## 2 28       under       7.45
## 3 35       canopy     10.1 
## 4 35       under       7.79